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ABSTRACT 


The  analysis  and  design  of  a  submarine  propulsor  requires  the  ability  to  pre¬ 
dict  the  characteristics  of  both  laminar  and  turbulent  flows  to  a  high  degree 
of  accuracy.  This  report  presents  results  of  certain  benchmark  computations 
based  on  an  upwind,  high-resolution,  finite- differencing  Navier-Stokes  solver. 

The  purpose  of  the  computations  is  to  evaluate  the  ability,  the  accuracy  and 
the  performance  of  the  solver  in  the  simulation  of  detailed  features  of  viscous 
flows.  Features  of  interest  include  flow  separation  and  reattachment,  sur¬ 
face  pressure  and  skin  friction  distributions.  Those  features  are  particularly 
relevant  to  the  propulsor  analysis.  Test  cases  with  a  wide  range  of  Reynolds 
numbers  are  selected;  therefore,  the  effects  of  the  convective  and  the  diffusive 
terms  of  the  solver  can  be  evaluated  separately.  Test  cases  include  flows  over 
bluff  bodies,  such  as  circular  cylinders  and  spheres,  at  various  low  Reynolds 
numbers,  flows  over  a  flat  plate  with  and  without  turbulence  effects,  and 
turbulent  flows  over  axisymmetric  bodies  with  and  without  propulsor  effects. 
Finally,  to  enhance  the  iterative  solution  procedure,  a  full  approximation 
scheme  V-cycle  multigrid  method  is  implemented.  Preliminary  results  indi¬ 
cate  that  the  method  significantly  reduces  the  computational  effort. 

ADMINISTRATIVE  INFORMATION 

This  investigation  was  authorized  by  the  Office  of  Naval  Research  under  6.2  Viscous 
Flow  Program,  in  accordance  with  Program  Element  62323N,  Task  Area  R2332MS3, 
and  Work  Request  Number  N001495WX20047/AJ.  This  work  was  performed  at  the 
Naval  Surface  Warfare  Center,  Carderock  Division  (NSWCCD),  David  Taylor  Model 
Basin  (DTMB)  under  Work  Unit  5060-567. 

INTRODUCTION 

Hydrodynamically  speaking,  propulsor  components  are  lifting  bodies  that  provide 
thrust  for  propulsion.  To  improve  the  propulsive  performance,  it  is  desirable  to  have  a 
lifting  body  with  an  optimum  lift  to  drag  ratio.  Successful  analysis  and  design  requires 
the  ability  to  predict  the  hydrodynamic  forces  on  the  lifting  body,  such  as  lift  and  drag, 
to  a  high  degree  of  accuracy.  Recently,  large  research  efforts  on  computational  fluid 
dynamics  (CFD)  have  been  directed  to  achieve  such  goals.  The  most  practical  approach 
is  to  derive  some  numerical  methods  for  the  Reynolds-averged  Navier-Stokes  equations 
(RANS).  The  more  popular  methods,  depending  on  the  manner  the  convective  terms 
are  formulated,  include  the  central  differencing  and  the  upwind  differencing  schemes. 
Both  upwind  (DTNS  Code1)  and  central  (IFLOW  Code2)  differencing  formulations  have 
been  used  for  hydrodynamics  simulations  at  the  David  Taylor  Model  Basin  (DTMB). 
As  part  of  DTMB  propulsor  research  efforts,  a  variation  of  the  upwind  schemes  has 
been  employed  to  simulate  the  flow  through  the  blade  rows  of  a  turbomachinery3,  4. 
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For  high  Reynolds  number  flows,  the  Navier-Stokes  equations  become  convectively 
dominated  and  are  hyperbolic  in  nature.  Under  such  circumstances,  the  upwind  differ¬ 
encing  approach  becomes  attractive  and  has  certain  advantages.  For  a  one-dimensional 
case,  it  can  be  shown  that  the  information  at  any  point  propagates  in  the  direction 
according  to  the  sign  of  the  local  eigenvalues  of  the  flux  Jacobian.  Consequently,  a 
numerical  procedure  with  the  upwind  differencing  technique  based  on  the  direction  of 
local  wave  propagation  is  physically  more  adaptive  to  the  characteristics  of  the  equa¬ 
tions.  Besides,  in  the  upwind  scheme,  the  dispersive  and  dissipative  errors  are  more 
closely  balanced  than  in  other  schemes  of  equal  or  even  higher  order  that  use  the  same 
set  of  nodal  points,  regardless  of  the  direction  of  the  convection5.  The  flux-differencing 
splitting  approach  suggested  by  Roe6  is  a  popular  upwind  scheme  for  solving  the  in¬ 
compressible  Navier-Stokes  equations,  since  it  does  not  require  the  inviscid  flux  to  be  a 
homogenous  function  of  order  one.  Based  on  Roe’s  approach6  the  flux-difference  is  split 
according  to  an  approximate  solution  of  the  Riemann  problem.  The  solution  provides 
the  information  about  the  direction  of  the  wave  propagation,  which,  in  turn,  is  incor¬ 
porated  into  the  discretized  system  to  form  an  upwind  scheme.  The  accuracy  of  the 
solution  can  be  promoted  to  a  higher  order  by  reconstructing  the  primitive  variables  or 
the  fluxes  midway  between  two  nodes  with  an  extrapolation  technique.  In  conjunction 
with  the  reconstruction  process,  the  slope  limiter  or  the  flux  limiter  can  be  implemented 
to  obtain  total  variation  diminishing  (TVD)  property7.  Viscous  terms  can  be  discretized 
with  a  second  order  central  differencing  scheme.  A  numerical  formulation,  based  on  the 
principles  described  above,  allows  the  discontinuity  of  the  solution  to  be  resolved  over 
only  two  adjacent  nodes  without  causing  the  non-physical  oscillations.  Therefore,  this 
approach  is  sometimes  called  high  resolution. 

This  report  presents  the  results  of  a  numerical  study  designed  to  evaluate  the  ability 
and  accuracy  of  an  upwind  scheme  in  predicting  certain  flow  features  that  are  relevant 
to  propulsor  analysis.  Some  of  the  features  of  interest  are  flow  separation  and  reattach¬ 
ment,  surface  pressure  and  skin  friction  distribution.  Test  cases  for  present  numerical 
study  are  selected  so  that  (1)  a  wide  range  of  Reynolds  numbers  are  covered,  (2)  the 
boundaries  are  smooth  and  the  distortions  of  the  grid  systems  are  minimized,  (3)  the 
flow  features  are  relevant  to  the  proplusor  analysis  and  (4)  the  test  data  are  well  ana¬ 
lyzed.  Test  cases  include  flows  over  bluff  bodies,  such  as  circular  cylinders  and  spheres, 
at  various  low  Reynolds  numbers,  flows  over  a  flat  plate  with  and  without  turbulence 
effects,  and  turbulent  flows  over  axisymmetric  bodies  with  and  without  propulsor  ef¬ 
fects.  To  enhance  the  iterative  solution  procedure,  a  full  approximation  scheme  (FAS) 
V-cycle  multigrid  method  is  implemented.  A  fast  convergence  rate  is  achieved  as  a 
result. 
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THE  GOVERNING  EQUATIONS 


The  three-dimensional  incompressible  RANS  equations  based  on  primitive  variables 
are  formulated  in  a  boundary-fitted  curvilinear  coordinate  system.  Using  Chorin’s  arti¬ 
ficial  compressibility  formulation,8  the  incompressible  Navier-Stokes  equation  is  written 
in  conservation  form  for  three-dimensional  flow  in  Cartesian  system  as 

Q,  +  (E'  -  El)r  +  (F*  -  F‘)„  +  (<?  -  G'J,  =  0  .  (1) 

In  Eq.  1,  the  dependent  variable  vector  Q  is  defined  as  Q  —  ( p ,  u,  v,  w)T ,  and  the  inviscid 
flux  vectors  E*,  F*,  and  G*,  the  viscous  shear  flux  vectors  E*,F*,  and  G*  are  given  by 

E*  =  (fiu,u2  +  p,uv,uw)T 
F*  =  (fiv,uv,v2  +  p,  vw)T 
G*  =  (f3w,uw,  vw,w2  +  p)T 
Ev  —  Rc  (0,  Uy  -f-  Uj,,  uz  -f  wx) 

F*  =  Re”1  (0,  uy  -f  vx,  2vy,  vz  +  wy)T 
G*  =  Re~1(0,uz  +  wx,vz  +  wy,2wz)T. 

The  coordinates  x,  y,  z  are  scaled  with  an  appropriate  characteristic  length  scale  L.  The 
Cartesian  velocity  components  u,v,w  are  nondimensionalized  with  respect  to  the  free 
stream  velocity,  Vo-  The  normalized  pressure  is  defined  as  p  =  (p  —  Poo)/pV£-  The 
kinematic  viscosity,  u,  is  assumed  to  be  constant,  and  the  Reynolds  number  is  defined  as 
Re  =  VooL/v.  The  artificial  compressibility  parameter,  /3,  monitors  the  error  associated 
with  the  addition  of  the  unsteady  pressure  term  dp/dt  in  the  continuity  equation.  The 
unsteady  pressure  term  is  needed  for  coupling  the  mass  and  momentum  equations  to 
make  the  system  hyperbolic. 

Equation  1  can  be  transferred  to  a  curvilinear,  body-fitted  coordinate  system  ((  ,  £  , 
rj  )  through  a  coordinate  transformation  of  the  form 

C  =  ({x,y,z),  £  =  ({x,y,z),  and  rj  =  y(x,y,z). 

Thus,  Eq.  1  becomes, 

(Q/J)t  +  (E-  Ev)<  +  (F-  Fv)t  +  {G-  G„)„  =  0  (2) 


with 

(E,F,G)T  =  {{T]  (E*/J  ,F*/J  ,G*/Jf} 

and 

(E„F„G,)T  =  [r]  (E-JJ  ,f;/j  ,crjjf  , 
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where 


m  = 


Cx  Cy  Cz 

£x  £y  £z 

L  Vx  %  rj2 


and  the  Jacobian  of  the  coordinate  transformation  is  given  by 


J  1  =  det 


x(  yc  zc 

xt  Vi  H 

xv  Vv  zv  . 


The  Jacobian  is  the  ratio  of  volume  elements  in  the  two  coordinate  systems.  For  a  proper 
transformation,  neither  J  nor  its  reciprocal  is  zero.  At  present,  the  transformation  is 
chosen  so  that  the  grid  spacing  in  the  computational  domain  is  uniform  and  of  one  unit 
in  length  in  all  three  spatial  directions.  The  Cartesian  derivatives  of  the  shear  fluxes  are 
obtained  by  expanding  them  using  chain  rule  expansions  in  the  (,  £,  and  r/  directions. 


NUMERICAL  DISCRETIZATION 

Let  i,j,k  denote  the  integer  indices  of  a  grid  point  in  the  curvilinear  system  (,£,r) 
whose  Cartesian  coordinates  are  x,y,z.  Each  grid  point  serves  also  as  the  centroid 
of  a  control  volume  whose  six  bounding  surfaces  are  formed  by  bisecting  the  distances 
between  the  centroid  and  its  six  adjacent  grid  points.  The  Cartesian  flow  variables  such 
as  u,v,w  and  p  are  placed  at  each  grid  point.  The  indices  of  the  grid  point  are  used 
as  subscripts  for  the  variables  to  indicate  the  association.  To  avoid  introducing  any  free 
stream  error,  the  metric  terms  such  as  (x,  £x  and  tjx  etc.,  are  computed  from  x,y,z 
data  by  using  a  second-order  central- differencing  approximation  of  x^,  x^  and  xv  etc., 
as  described  by  Pulliam  and  Steger9. 

INVISCID  FLUXES 

At  present,  discretization  of  the  inviscid  fluxes  of  Eq.  2  is  achieved  by  applying 
the  Riemann  solver  to  each  direction  of  the  coordinate  system  as  suggested  by  Yee, 
Warming  and  Harten10.  Consider  a  one-dimensional  hyperbolic  system  of  conservative 
laws 


(§),+ *«)“«■  =  w  (3) 

where  6  =  (,  £  or  T),  and  D(Q )  is  the  Jacobian  Matrix.  Both  Q  and  H(Q)  are  column 
vectors  with  four  components  and  D(Q )  is  a  four  by  four  matrix. 

Let  the  right  eigenvectors  be  the  column  elements  of  matrix  R  and  eigenvalues 
Ai,  A2,  A3  and  A4  be  the  diagonal  elements  of  matrix  A,  then  the  relationship  between 
D,  R  and  A  is  given  by  a  similarity  transformation  D  =  RAR*1 .  The  row  elements 
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of  the  matrix  give  an  orthonormal  set  of  the  left  eigenvectors.  In  the  discretized 
system,  the  state  Qi  at  grid  point  l  is  considered  as  an  averaged  value  in  an  interval, 
that  is 

/•(Z+i)  A0 

Qi  =  A r1  /  Q  dO.  (4) 

Roe’s  flux  difference  splitting  is  constructed  by  forming  a  mean  value  Jacobian  D(Qi+ 1,  Qi) 
such  that 

H(Qi+ 1)  -  H(Q,)  =  D(QI+UQ,)  (qw  -  (?,), 

and 


D(Qi+i,  Qi)  —  D(Qi+i,  Qi)  as  Qi+i Qi  ■ 

to  evaluate  the  mean  value  (the  local  frozen  value  )  Jacobian  at  the  interface  /  ±  |,  the 
simple  average  value  of  the  state  Q  at  two  adjacent  grid  points  are  used,  that  is 

A±!  =  A±§(Qz±§)  (5) 

where 

Qi±±  =  t^Qi  ±  <5z±i) 

By  the  similarity  transformation,  Dl+ 1  can  be  written  as 

A+|  =  (RkR~i)i+i 

=  (RA+ R_1)l+i  —  (RA~  R~')l+i 
=  A+i  ~  &1- 1  ? 

where  A+  and  A-  are  the  absolute  values  of  the  positive  (speed  of  the  right  travelling 
wave)  and  the  negative  (speed  of  the  left  travelling  wave)  eigenvalues,  respectively. 
With  the  mean  value  Jacobian  locally  defined,  Eq.  3  can  be  decoupled  into  a  system  of 
four  scalar  equations  with  the  eigenvalues  representing  the  wave  speeds  of  the  Riemann 
problem.  Let’s  define  the  local  characteristic  variables  W  as  the  projection  of  Q  into 
the  left  eigenvectors  R-1  ;  therefore  ,  W  —  R~'Q.  For  any  given  two  states  Qi  and 
Qi+i,  the  flux  at  the  interface  Hl+ r  can  be  expressed,  in  term  of  flux  difference,  as  either 

Hl+>_  =  Ht  +  RA~  Al+,W  (6) 

or 

H,H=Hw-RA*A,+iW  ,  (7) 
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where  Al+±W  =  (WJ+1  -  Wi)/2  . 

Assuming  that  R ,  R  1  and  A  are  constant,  with  Eqs  6  and  7,  Eq.  3  can  be  written  as 
a  system  of  four  scalar  hyperbolic  conservative  equations  for  the  characteristic  variables 
W  ,  that  is, 


fWm\ 

\-J- J  ~  +  (A+),_iA =  0  for  m  =  1,2, 3, 4. 


An  implicit  delta  form  of  Eq.  8  is 


1 


tJ  ~  (Am)/+iA/+i  +  (A+);_iA,_i 


A  W” 


=  (A^WA^WS  -  (A+)f_xAMW^  , 


(8) 


(9) 


where  r  is  the  time  step  size,  n  is  the  time  step  number,  and  AW”  =  W”+1  -  W”. 

When  W”  is  sufficiently  differentiable,  the  local  jump  in  W”  at  the  interfaces  l  ±  |, 
known  as 


A, ,  i  Wn 

(±5  TO 

(10) 

are  replaced  with 

(11) 

Equation  10  is  only  a  first-order  one-sided  approximation  of  Eq.  11.  To  enhance  the 
accuracy,  the  following  relationships  from  a  scalar  scheme  can  be  extended  to  a  constant 
coefficient  system  by  applying  them  scalarly  to  each  of  the  m  scalar  characteristic 
components  in  Eq.  8, 


and 


'3Wm\” 


V  36 


_A6  =  (w£,,  - 


+ 


I[(l  +w)(4-,I+1  -  4>t,,)Al+iW^  +  (1  ) A, -lW: 


,(12) 


-j[(l  +“)(*;,/  -  +  (1  -u)(Kj»  -  tJ.,)A,+JW“]  ,(13) 

for  m=l,2,3,4,  with 


(14) 
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and 


±  f  (AH^"/AI+iI^)±  forA,±jH^O 

l»  forA,±iH"=0  '  (15) 

The  order  of  the  accuracy  in  the  spatial  derivatives  presented  above  is  determined  by  the 
values  of  u.  For  u  =  —  1,  the  scheme  is  fully  upwind  second-order  accurate.  For 
the  scheme  is  upwind  biased  third-order  accurate.  Function  $  is  called  the  limiter;  it 
is  used  to  control  unwanted  oscillations  in  numerical  schemes.  Various  designs  of  the 
limiter  were  found  and  successfully  tested.  3’ 11 

Substitute  Eqs.  12  and  13  into  Eq.  9,  and  multiply  Eq.  9  by  the  set  of  right  eigen¬ 
vectors  R  from  the  left.  A  conservative  high-resolution  scheme  for  the  nonlinear  system 
is  derived: 

with 


Ar±i  -  io  - ")Aai + o + “)Af±i](»<±i  -  )/4l” , 

.  2  2  2 

and 


-  (D-.A^-D+A^r 


(RK-'R~%,  Al+hQn 


vz+§' 


A  Qn  = 

(RK+1  R-1)?_iAl_iQn 


(16) 


The  right  hand  side  of  the  Eq.  16  is  evaluated  at  time  level  n;  it  is  the  spatial 
derivative  of  Eq.  2  and  is  designated  as  residual.  Equation  16  represents  the  relationship 
between  the  residual  at  nth  time  step  and  the  correction  of  the  solutions  at  n  -f  lt/l  time 
step.  The  correction  and  the  residual  approach  to  zero  as  the  solutions  approach  to 
their  steady  state  values. 

VISCOUS  FLUXES 

The  viscous  fluxes  in  Eq.  2  are  evaluted  by  a  second-order  central  differencing 
scheme.  The  computation  of  the  fluxes  require  all  nine  metric  coefficients  at  each 
of  the  six  bounding  surfaces  of  each  computational  cell. 

SOLUTION  ALGORITHM 

Equation  16  can  be  extended  to  three-dimensional  applications  with  the  operator 
split  method.  The  differencing  schemes  described  previously  are  applied  to  each  co¬ 
ordinate  direction  (,  £,  and  rj  independently;  a  summation  over  all  directions  gives 
the  discretized  approximation  of  a  multi-dimensional  flow  problem.  Upon  forming  the 
Jacobian  matrices  A ,  B ,  and  C  from  invicid  fluxes  E,  F,  and  G  and  X ,  Y ,  and 
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Z  from  the  viscous  fluxes  Ev,  Fv,  and  Gv  ,  a  finite- differencing  form  of  Navier-Stokes 
equations  can  then  be  written  as: 

(i_  -  (A-+X)i+JAj+i  +  (A+  +  X)j.iA,., 

-  (B-  +  y)J+iAj+i  +  (B*  +  nH4H 

-  (C-  +  Z)i+i\H  +  (C-  +  AQ-  =  -RESm  ■  (U) 

Equation  17  is  solved  by  an  implicit  hybrid  algorithm  where  a  symmetric  planar 
Gauss-Seidel  relaxation  is  used  in  the  (  direction  and  approximation  factorization  is 
applied  in  the  £  and  77  directions: 


M  -  (B~  +  K)i+j Ai+i  +  (B+  +  n- jAi-i 
M  —  (C  +  Z)k+iAk+i  +  (C+ +  Z)k_iAk_i 


AQ  =  -RES(Qn,Qn+1)  , 
AQ  =  MAQ  , 


(18) 

(19) 


gn+1  =  +  AQ"  ,  (20) 

with  M  =  //(rJ)  +  ( A~  +  X)i+i  +  (A+  +  X)-_i,  where  RES(Qn,  <2n+1)  indicates  the 
nonlinear  updating  of  the  residual  while  sweeping  in  the  £  direction. 


MULTIGRID  METHOD 

For  certain  simulations,  in  order  to  obtain  meaningful  results,  a  large  number  of 
finely  sized  grid  points  is  needed.  The  adverse  effect  of  such  a  grid  system  upon  the 
computation  effort  is  that  the  rate  of  convergence  deteriorates  significantly.  The  analysis 
provided  by  Brandt12  suggests  that  the  low-frequency  components  of  the  errors  are 
efficiently  approximated  on  coarse  grids  but  are  slow  to  convergence  on  fine  grids. 
In  addition,  the  high-frequency  components  must  be  approximated  on  fine  grids.  By 
utilizing  interactively  several  scales  of  discretization,  multigrid  techniques  resolve  such 
conflicts  and  avoid  stalling. 

To  accommodate  nonlinearities,  a  full  approximate  scheme  (FAS)  is  used.  The 
discretized  system  of  equations  described  previously  can  be  represented  as: 


HQ)  =  -R , 


(21) 


where  L  is  the  differencing  operator,  Q  is  the  unknown  to  be  solved  and  R  is  the 
residual.  The  iterative  process  will  reduce  the  residual  to  zero  as  the  steady-state 
solution  is  approached.  The  FAS  procedures  for  solving  Eq.  21  can  be  described  as 
follows: 

(1)  relax  on  the  fine  grid,  Lh{Qh)  =  —Rh  , 
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(2)  solve  L2h{Q2h)  +  [Ilh[Lh{Qh)}  -  L2h(I2hQh)  =  ~R2h  on  the  coarse  grid,  and 

(3)  replace  Qh  < —  Qh  +  I%h(Q2h  —  I2hQh)  on  the  fine  grid. 

The  notation  introduced  here  includes  the  difference  operators  at  the  fine  grid  Lh 
and  the  coarse  grid  L2h  ,  the  restriction  operators  I2h  (for  the  approximation)  and 
I2h  (for  the  residual),  and  the  interpolation  operator  I^h. 

RESULTS 

In  the  followings  benchmark  computations,  in  order  to  qualify  the  comparison  be¬ 
tween  the  experimental  measurements  and  computational  results,  root-mean-square 
(RMS)  differences2  are  calculated.  The  RMS  difference  is  defined  as 

RMS  =  [XK  -  oD’/Jv] 5  . 

L  i  J 

where  vc  and  vm  are  computed  and  measured  values  respectively,  N  is  the  total 
number  of  data  values  used  in  the  comparison,  and  the  subscript  i  ranges  from  1  to 
N.  Computations  were  carried  out  with  64-bit  precision  on  a  Silicon  Graphics  Power 
Challenge  machine. 

TURBULENT  FLOW  ALONG  A  FLAT  PLATE 

The  structure  of  the  turbulent  boundary  layer  along  a  flat  plate  has  been  investigated 
earlier  by  Ludweig  and  Tillmann13,  and  Wieghart  and  Tillmann14.  It  was  found  that 
the  axial  velocity  profile  in  the  inner  one-fifth  of  the  boundary  layer  can  be  represented 
by  the  universal  logarithmic  law  (excluding  the  viscous  sublayer).  The  remaining  outer 
four-fifths  can  be  adequatedly  expressed  by  the  power  law.  Wieghart  and  Tillmann’s 
data  were  collected  in  a  wind  tunnel  test.  The  flow  velocity  was  33  m/sec  and  the 
average  dynamic  viscosity  was  0.151  cm2/sec.  Velocity  measurements  were  taken  at 
twenty-three  locations  ranging  from  0.087  m  to  4.987  m  from  the  leading  edge.  The 
boundary  thickness  grew  from  0.0335  cm  to  0.9242  cm.  The  test  data  were  compiled 
and  presented  at  a  1968  Stanford  turbulent  flow  conference15. 

In  the  present  computation,  the  computational  domain  extends  8  m  in  streamwise 
direction,  0.16  m  in  cross  flow  direction,  and  0.5  m  in  the  third  direction;  grid  points 
used  are  57,  61  and  5  in  the  respective  directions.  The  Reynolds  number  (Re)  is  2.2x1 06. 
Grid  distribution  in  the  cross-flow  direction  is  non-uniform,  and  is  clustered  near  the 
plate  such  that  the  y+  coordinate  of  the  first  grid  point  off  the  plate  surface  is  less 
than  0.3.  The  Courant-Friedrichs-Lewy  (CFL)  number  for  present  computations  is 
103.  For  the  turbulence  modelling,  the  standard  Baldwin  and  Lomax’s  algebraic  eddy 
viscosity  formulation  16  is  used.  Figure  la  shows  the  skin  friction  coefficient  Cr  along 
the  plate  surface,  and  Figs,  lb  and  lc  show  the  axial  velocity  profiles  at  x=0.78  m  and 
4.98  m,  respectively.  The  RMS  differences  indicate  that  the  deviations  between  the 
measurements  and  the  predictions  are  within  the  limits  of  the  expected  measurement 
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uncertainties.  Figure  2  shows  the  velocity  profile  at  x=4.98  m  presented  in  (y+ ,  u+) 
coordinates.  It  can  be  seen  that  the  predictions  agree  well  with  the  measurements.  In 
the  turbulent  zone  where  y+  lies  between  30  and  1000,  both  the  predictions  and  the 
measurements  fit  the  universal  logarithmic  law.  However,  none  of  the  data  lie  within 
the  sublayer.  The  convergence  history  presented  in  Fig.  3  shows  that  the  residual  value 
approaches  the  machine  zero. 


X(m) 


Fig.  la.  Skin  friction  coefficient. 


Fig.  lb.  Velocity  profile  at  x=0.78  m.  Fig.  lc.  Velocity  profile  at  x=4.98  m. 

Fig.  1.  Skin  friction  coefficient  and  velocity  profile 
for  turbulent  flow  along  a  flat  plate. 


10 


Fig.  2.  Velocity  profile  in  (y+,u+)  coordinates. 


Fig.  3.  Convergence  histories  for  computing 
turbulent  flow  along  a  flat  plate. 


LAMINAR  FLOW  ALONG  A  FLAT  PLATE 


For  steady  laminar  flow  along  a  flat  plate  with  a  constant  free  stream  velocity 
U0  ,  the  pressure  gradient  dp/ dx  in  streamwise  direction  vanishes.  The  Navier-Stokes 
equations  reduce  to  the  Prandtl  boundary-layer  equations.  A  solution,  known  as  the 
Blasius  solution  after  its  originator,  is  obtained  by  assuming  similar  profiles  along  the 
plate  at  every  location  along  the  plate.  Blasius  assumed  that 


where  y  is  the  distance  above  the  plate  surface,  8  is  the  boundary  layer  thickness,  and 


y  y 

6  ~  x/s/R,  ~  V  ’ 


with  x  the  distance  from  leading  edge  and  Rx  the  Reynolds  number  based  on  x.  Under 
the  similarity  assumption,  the  Prandtl  boundary  layer  equations  can  be  further  reduced 
to  an  ordinary  differential  equation.  Solution  can  then  be  obtained  numerically.  In  the 
present  numerical  simulation,  the  geometrical  dimension  used  previously  for  turbulent 
flow  is  adopted  and  the  Reynolds  number  is  3.64xl05.  The  grid  size  is  129x129x3. 
Figure  4a  shows  the  skin  friction  coefficient  along  the  plate  surface,  and  Figs.  4b  and 
4c  show  the  profiles  of  axial  velocity  u/U0  ,  and  transverse  velocity  w/U0  ,  respectively. 
The  RMS  difference  for  each  quantity  is  also  presented.  The  results  of  the  present 
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Fig.  4c.  Transverse  velocity  profile. 


Fig.  4.  Skin  friction  coefficient  and  velocity  profiles 
for  laminar  flow  along  a  flat  plate. 


computation  agree  well  with  the  solutions  obtained  by  Blasius’s  method.  Part  of  the 
deviation  can  be  attributed  to  the  fact  that  the  Blasius’s  solution  is  based  on  Prandtl 
boundary  layer  equations  while  the  current  computation  is  based  on  Navier-Stokes 
equations.  Figure  5  shows  the  convergence  histories  for  solutions  with  multigrid  (7 
levels)  and  without  multigrid  (1  level)  application.  Considerable  computing  effort  is 
saved  with  the  application  of  the  multigrid  method. 
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Fig.  5.  Convergence  histories  for  computing 
laminar  flow  along  a  flat  plate. 


FLOW  PAST  CIRCULAR  CYLINDERS  AT  LOW  SPEEDS 

For  flow  past  a  circular  cylinder  at  Reynolds  numbers  (length  scale  is  based  on  the 
diameter)  below  the  critical  value  (  ~  40)  at  which  a  Karman  street  is  formed,  the  flow 
is  steady  and  twin  vortices  exist  behind  the  cylinder.  There  are  many  interesting  flow 
related  phenomena  despite  the  simplicity  of  the  geometry  involved.  The  phenomena 
are  the  boundary  layer  development  along  curved  surface,  the  flow  separation,  and 
wake  reattachment.  However,  the  details  of  these  phenomena,  such  as  the  locations  of 
the  separation,  the  coordinates  of  the  vortex  center,  and  the  wake’s  length  and  shape, 
are  Reynolds  numbers  dependent  and  pose  great  challenges  for  numerical  simulation! 
For  these  and  other  reasons,  this  problem  has  drawn  much  attention  in  the  past,  both 
theoretically  and  experimentally.  In  light  of  the  availability  and  quality  of  the  data, 
this  problem  is  selected  as  a  benchmark  case  for  present  numerical  study. 

The  outer  boundary  of  the  computational  domain  is  described  by  a  circle  with  a 
radius  10  times  that  of  the  cylinder.  An  O-type  and  orthogonal  grid  system  is  se¬ 
lected.  The  grid  distribution  in  the  radial  direction  is  non-uniform  and  is  clustered 
near  the  surface.  The  distance  between  the  cylinder  surface  and  the  first  grid  point  is 
one-thousandth  of  the  radius  of  the  cylinder.  For  computation,  symmetry  is  assumed 
and  only  the  plane  above  the  axis  of  symmetry  is  considered.  The  number  of  grid 
points  in  radial,  circumferential  and  axial  directions  are  257,  257  and  3,  respectively. 
The  CFL  number  used  for  the  following  computations  is  102.  The  boundary  conditions 
are:  (1)  a  non-slip  ,  non-penetrating,  and  vanishing  normal  pressure  gradient  on  solid 
surface,  (2)  prescribed  free  stream  values  at  the  upstream  side  of  the  outer  boundary, 
(3)  second-order  extrapolations  at  the  downstream  side  of  the  outer  boundary  and  (4) 
periodicity  in  the  spanwise  direction.  In  order  to  describe  the  main  features  of  the  flow, 
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Fig.  6.  Geometrical  parameters  of  the  closed  wake. 


Coutanceau  and  Bernard’s17  geometrical  parameters  for  the  closed  wake  are  adopted 
and  shown  in  Fig.  6,  where  parameters  l  and  L  are  the  width  and  length  of  the  wake,  6S 
is  the  separation  angle  and  the  vortex  centers  are  located  on  the  ( x,y )  axes  by  a  and 
b  .  It  was  observed  by  Taneda18  that  the  twin  vortices  behind  the  cylinder  appeared 
when  the  Reynolds  number  was  greater  than  6  and  became  unstable  when  the  Reynolds 
number  was  greater  than  45.  Therefore,  the  Reynolds  numbers  5  and  40  are  selected  as 
the  lower  and  upper  bounds  for  the  present  steady-state  computations.  Taneda’s18  and 
Coutanceau  and  Bouard’s17  flow  visualizations  were  obtained  by  similar  methods:  by 
illuminating  the  light  particles  suspended  uniformly  in  the  liquid  and  by  photographing 
in  the  direction  normal  to  the  lighted  plane.  Taneda18  used  aluminium  and  water  in 
the  tests,  while  Coutanceau  and  Bouard17  used  fine  bright  particles  and  the  Vaseline 
oil  ‘  MARCOL  80  ’  in  their  tests.  The  latter  reference  derived  the  particle  velocity 
by  measuring  the  length  of  the  particle  trajectory  during  the  time  of  exposure.  The 
reported  inaccuracy  was  less  than  2%.  The  wall  influence  was  investigated  by  chang¬ 
ing  the  ratio  A  between  the  cylinder  and  the  tank  diameters17.  Figure  7  shows  the 
computed  wake  shapes  behind  the  cylinder  at  different  Reynolds  numbers.  Compared 
with  the  flow  visualizations17, 18 ,  the  characteristics  of  the  wake  shapes  near  the  sepa¬ 
ration  and  reattachment  points  are  well  simulated.  The  result  indicates  that  the  twin 
vortices  begin  to  develop  at  Reynolds  number  about  5  and  it  agrees  with  Taneda’s18 
observation.  In  Fig.  8,  the  currently  computed  separation  angles  at  various  Reynolds 
numbers  are  compared  with  those  computed  earlier  by  Keller19’  20  and  Raal21  and  those 
measured  by  Coutanceau  and  Bouard.17  Figure  9  shows  the  relationship  between  the 
Reynolds  numbers  Re  and  the  wake  lengths  L/D,  from  both  the  current  computations 
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Fig.  8.  Separation  angles  at  different  Reynolds  Fig.  9.  Wake  lengths  at  different  Reynolds 
numbers  for  flow  past  a  circular  cylinder.  numbers  for  flow  past  a  circular  cylinder. 

and  Taneda’s  observations.18  In  Fig.  10,  the  coordinates  of  the  vortex  center  (a,b)  (see 
Fig.  6)  are  plotted  against  the  Reynolds  number,  Re.  The  experimental  data  were  re¬ 
ported  by  Coutanceau  and  Bouard17  with  A  ,  the  ratio  between  the  cylinder  and  the 
tank  diameters,  equal  to  0.024.  The  small  value  of  A  indicates  that  the  wall  effect  is 
relatively  small.  Figure  11  shows  the  similarity  of  the  closed- wake  shape.17  For  com¬ 
puted  and  measured  wakes  at  different  Reynolds  numbers,  when  the  wake  width  /  and 
its  distance  from  the  rear  stagnation  point  X  —  R,  are  normalized  with  the  maximun 
wake  width  lmax  and  wake  length  L,  respectively,  and  then  plotted  against  each  other, 
the  results  merge  into  a  single  curve  except  at  the  regions  near  the  cylinder  wall.  Fig¬ 
ure  12  shows  the  velocity  similarity  on  the  rear  flow  axis  in  the  closed  wake,  where  the 
velocity  u  is  normalized  with  its  maximum  value  umax.  Figure  13  shows  the  compar¬ 
isons  of  the  computed  and  measured  velocities  at  rear  flow  axis  at  Reynolds  numbers 
20  and  40.  The  RMS  differences  are  comparable  with  the  measurement  uncertainty 
(~  2%).  Figure  14  shows  the  computed  pressure  distribution  at  the  cylinder  surface 
for  Re=40;  it  is  compared  with  the  distributions  observed  at  Re=36  and  Re=45  by 
Thom,22  and  computed  by  Apelt23  at  Re=40.  Thom’s  approximate  theory24  for  deter¬ 
mining  the  value  of  the  pressure  at  the  front  stagnation  point  at  low  speed  gives  the 
result  at  Reynolds  number  40,  (1+7/Re)  or  1.175,  which  agrees  well  with  the  value  of 
the  current  result  1.18.  Figure  15  presents  the  convergence  histories  of  the  numerical 
simulation  at  Re=40.  The  7-level  multigrid  solver  improves  the  efficiency  significantly. 
A  fine  grid  size  is  needed  for  detailed  computation,  because  with  such  a  fine  grid  size 
the  convergence  is  slow.  The  application  of  multigrid  technique  (7  levels)  reduces  the 
computational  effort  considerably  for  a  given  CPU  time. 
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a/D,  1/2  x  b/D 


Reynolds  number.  Re  (X-R)/L 


Fig.  10.  Locations  of  vortex  centers  of  wakes  Fig.  1 1.  Similarity  of  the  closed-wake  shapes, 
at  different  Reynolds  numbers. 


(X-R)/L 


Fig.  12.  Velocity  similarity  on  the  rear  flow  axes 
of  the  closed  wakes. 
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fl/n 


Fig.  13.  Velocity  distributions  on  the  flow  axis  behind  the  circular  cylinder. 


Fig.  14.  Surface  pressure  distribution 
for  flow  past  a  circular  cylinder  at  Re=40. 


Fig.  15.  Convergence  histories  for  computing 
flow  past  a  circular  cylinder  at  Re=40. 
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FLOW  PAST  SPHERES  AT  LOW  SPEEDS 


The  flow  and  the  solutions  are  assumed  to  be  axisymmetric.  The  topology  of  the 
grid  and  the  specification  of  the  boundary  conditions  are  the  same  as  those  used  for 
computing  the  flow  over  the  cylinder  presented  earlier,  except  that  the  reflective  con¬ 
dition  is  applied  in  the  circumferential  direction.  The  solutions  on  different  meridional 
planes  are  related  by  simple  coordinate  transformation.  The  CFL  number  used  for  the 
computations  is  102.  Figure  16  shows  the  computed  wake  shapes  behind  the  sphere  at 
various  Reynolds  numbers.  The  relationship  between  the  Reynolds  number  Re  and  the 


Fig.  16.  Standing  vortices  downstream  of  a  sphere. 
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wake  length  L  is  plotted  in  Fig.  17;  also  included  are  the  test  data  from  Taneda.25 
The  twin  vortices  behind  the  sphere  appear  at  Reynolds  number  about  25.  The  surface 
pressure  distribution  at  Reynolds  number  100  is  shown  in  Fig.  18;  also  included  are  the 
results  from  earlier  computations.26,  27  The  convergence  histories  for  Re=100  is  shown 
in  Fig.  19.  Multigrid  technique  improves  the  computational  efficiency  significantly. 


Reynolds  number.  Re 

Fig.  17.  Wake  lengths  at  different  Reynolds 
numbers  for  flow  past  a  sphere. 


0  50  100  150 

Degree  from  front  stagnation  point 


Fig.  18.  Surface  pressure  distribution 
for  flow  past  a  sphere  at  Re=100. 


2000  4000  6000 

Woik  unit 

Fig.  19.  Convergence  histories  for  computing 
flow  past  a  sphere  at  Re=100. 
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TURBULENT  FLOW  PAST  AXISYMMETRIC  BODIES 

The  purpose  of  this  computation  is  to  investigate  the  accuracy  of  the  numerical 
scheme  in  predicting  the  characteristics  of  the  flow  in  an  expanding  boundary  layer 
under  adverse  pressure  gradient  at  high  Reynolds  numbers.  Axisymmetric  bodies,  des¬ 
ignated  as  DTMB-bodyl  and  DTMB-body2,  were  built  and  tested  at  the  David  Taylor 
Model  Basm  (DTMB)28.  The  Reynolds  number  at  the  test  condition  was  6.6xl06.  The 
size  of  the  C-type  grid  for  the  present  computation  is  197x3x146  in  radial,  circumfer¬ 
ential  and  axial  directions,  respectively.  Modified  Baldwin-Lomax  model2  is  used  to 
simulate  the  turbulent  flow.  Figure  20  shows  the  computed  (solid  line)  and  measured 
(symbol)  pressure  distributions  on  surface  of  DTMB-bodyl.  The  RMS  difference  is 
0.019.  The  measurement  uncertainty  for  pressure  was  ±0.015.  It  can  be  seen  that 
the  pressure  gradient  is  zero  at  the  parallel  mid-body  region,  and  adverse  gradient  is 
present  at  stern  region.  In  Fig.  21  the  computed  velocity  profiles  at  various  axial  loca¬ 
tions  at  stern  region  of  DTMB-bodyl  are  compared  with  the  measurments.  The  RMS 
differences  range  from  0.012  to  0.025.  The  measurement  uncertainty  for  velocity  was 
±0.025.  In  Fig.  22  the  velocity  profiles,  at  the  parallel  mid-body  section,  are  plotted  in 
the  (y+,  u+)  coordinates.  The  value  of  the  y+  of  the  grid  nearest  to  the  body  surface 
is  about  2.5.  There  are  four  grid  points  that  lie  within  the  laminar  sublayer.  The 
computed  distribution  of  skin  friction  (solid  line)  on  body  surface  is  shown  in  Fig.  23. 
Compared  with  the  measurements,  the  RMS  difference  is  0.00029.  The  measurement 
uncertainty  for  skin  friction  was  ±0.0002.  The  computed  (solid  line)  and  measured 
(symbol)  turbulence  shear  stresses  near  the  stern  region  at  several  axial  locations  are 
shown  in  Fig.  24;  the  RMS  differences  range  from  0.008  to  0.015.  The  measurement 
uncertainty  for  turbulent  shear  stresses  is  ±0.01. 


Fig.  20.  Surface  pressure  distribution 
for  DTMB-bodyl. 
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Fig-21.  Velocity  profiles  at  stem  region  of  DTMB-bodyl 
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Fig.  22a.  x/L=0.4473 

Fig.  22b.  x/L=0.5765 

Fig.  22.  Comparisons  of  wall  law  at  mid-body  section  of  DTMB-bodyl 
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Fig.  23.  Skin  friction  coefficient 
along  the  surface  of  DTMB-bodyl 


Fig.  24.  Turbulent  shear  stresses  near  the  stem  region  of  DTMB-bodyl. 
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A  propeller  was  placed  on  DTMB-bodyl  whose  center  line  is  located  at  x/L=0.983. 
The  ratio  of  propeller  and  body  diameter  is  0.54.  Velocity  measurements  were  taken 
at  a  distance  of  0.227  propeller  diameter  upstream  of  the  propeller.  For  numerical 
prediction,  the  propeller  effect  was  simulated  with  a  body  force  model.  The  computed 

and  measured  velocity  profiles  under  two  different  propeller  operating  conditions  are 
shown  in  Fig.  25. 

The  same  type  of  computations  were  carried  out  on  DTMB-body2.  The  velocity 
profiles  at  different  axial  locations  are  shown  in  Fig.  26.  The  RMS  differences  between 
the  computed  and  measured  values  range  from  0.008  to  0.058. 


Fig.  25.  Velocity  profiles  upstream  of  an  operating  propeller 
at  different  advance  ratios. 
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Fie.  26a.  x/L=0.840. 


Fig.  26.  Velocity  profiles  at  stem  region  of  DTMB-body2. 
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CONCLUSIONS  AND  RECOMMENDATIONS 

An  upwind,  high-resolution,  finite- differencing  Navier-Stokes  solver  is  used  to  sim¬ 
ulate  steady-state  incompressible  flows  with  a  wide  range  of  Reynolds  numbers.  Ex¬ 
periments  with  measurements  of  high  quality  are  selected  as  benchmark  cases.  The 
predictions  are  compared  with  the  measurements  by  evaluating  the  root  mean  square 
(RMS)  differences.  In  all  cases,  the  RMS  differences  are  compatible  with  the  mea¬ 
surement  uncertainties.  For  the  low  Reynolds  number  cases,  the  detailed  features  of 
the  standing  vortices  behind  the  bluff  bodies  are  successfully  simulated.  For  the  high 
Reynolds  number  cases,  the  skin  friction  coefficients,  the  structure  of  the  turbulent 
velocity  profiles  and  the  turbulent  shear  stresses  are  correctly  predicted. 

At  low  Reynolds  numbers,  the  flows  are  dominated  by  diffusive  process.  The  rate  of 
convergence  of  the  iterative  procedure  becomes  very  slow,  even  for  an  implicit  method 
with  a  rather  high  CFL  number.  The  situations  can  be  improved  significantly  by 
implementing  the  multigrid  technique.  Compared  with  the  single-grid  approach,  the 
multigrid  solver  is  rather  insensitive  to  the  CFL  numbers  and  an  order  of  magnitude  of 
Central  Processor  Unit  (CPU)  time  is  saved. 

In  summary,  the  flow  features  that  are  relevant  to  the  propulsor  analysis,  such  as 
flow  separation  and  reattachment,  surface  pressure  and  skin  friction  dstributions,  can 
be  predicted  accurately  and  efficiently  with  an  upwind  RANS  solver.  The  formulations 
of  DTNS1  are  similar  to  the  formulations  described  in  this  report.  It  is  expected  that 
DTNS1  code  may  achieve  the  similar  performaces.  For  complicated  turbulent  flows,  to 
correctly  predict  the  turbulent  structures,  a  sophisticated  non-equilibrium  turbulence 
model  is  needed. 
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